Fluctuations and Correlations in Lattice Models for Predator-Prey Interaction 
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Including spatial structure and stochastic noise invalidates the classical Lotka-Volterra picture of 
stable regular population cycles emerging in models for predator-prey interactions. Growth-limiting 
terms for the prey induce a continuous extinction threshold for the predator population whose critical 
properties are in the directed percolation universality class. Here, we discuss the robustness of this 
scenario by considering an ecologically inspired stochastic lattice predator-prey model variant where 
the predation process includes next-nearest-neighbor interactions. We find that the corresponding 
stochastic model reproduces the above scenario in dimensions 1 < d < 4, in contrast with mean-field 
theory which predicts a first-order phase transition. However, the mean-field features are recovered 
upon allowing for nearest-neighbor particle exchange processes, provided these are sufficiently fast. 

PACS numbers: 87.23.Cc, 02.50.Ey, 05.40.-a, 05.70.Fh 



In 1920 and 1926, respectively, Lotka [l| and Volterra 
Q devised a simple coupled set of differential equa- 
tions to describe an autocatalytic reaction model and 
the statistics of fish catches in the Adriatic. The Lotka- 
Volterra model (LVM) has since become one of the cen- 
tral paradigms for the emergence of periodic oscillations 
in nonlinear systems with competing constituents Q , and 
features prominently in textbooks, from undergraduate- 
level population biology [i| to ecology ya, Ig and math- 
ematical biology as, for instance, it can also be for- 
mulated as a host-pathogen model j8j- Yet it has often 
been severely criticized as beingbiologically unrealistic 
and mathematically unstable pTul 13 • Recent investiga- 
tions of zero-dimensional 0] an d spatial stochastic mod- 
els H [H [H [H Q E3 have shown that this criticism 
definitely pertains to the original deterministic rate equa- 
tions; however, it turns out that the stochastic, or lattice, 
two-species predator-prey model variants display quite 
robust properties, rather insensitive on the details of the 
underlying microscopic processes (for a recent overview, 
see Ref. |l6|). In particular, the lattice predator-prey 
models (LPPM), display the following features: (i) The 
population densities typically display erratic (rather than 
regular periodic) oscillations, with amplitudes that van- 
ish in the thermodynamic limit ^| , caused by persistent 
and recurrent predator-prey activity waves that form 
complex spatio-temporal structures [l7|: (ii) when the 
prey population growth is limited (finite carrying ca- 
pacity, local site restrictions), there exists an extinction 
threshold for the predator population [l4l Il5| ; this con- 
stitutes a noncquilibrium active-to-absorbing-state phase 
transition with the critical exponents of directed percola- 
tion (DP) [l8L ll!^|. Also, for host-pathogen models with 
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two types of pathogens, the invasion of the system by one 
pathogen (the other becoming extinct) through oscilla- 
tory behavior, was reported using mean-field and pair- 
approximation treatments Q . 

As noted by various authors 0,0,0,0], a more re- 
alistic description of the predator-prey interaction should 
include the possibility for the agents to move. In fact, 
in real ecosystems prey tend to evade the predators, 
while the predators aim to pursue the prey. One ap- 
proach to account for the motion of the agents is to con- 
sider diffusion (nearest-neighbor hopping) of predators 
and/or prey, which however does not really affect the 
global properties of the LPPM [l4j . Another approach, 
to be considered here, is to assume a nearest-neighbor ex- 
change process (among any two agents: predators, prey 
and empty sites) in the following referred to as 'stir- 
ring'. It has to be noticed that both diffusion and stir- 
ring processes are not taken into account at the (mean- 
field) rate equation level. In addition, some recent in- 
vestigations have included long-range processes in two- 
dimensional LPPM, reporting quite different results on 
the existence 0] or absence 0] of 'self-sustained oscilla- 
tions' (in the thermodynamic limit). We also notice that 
for spatial host-pathogen models (with nearest-neighbor 
interactions), the rate equations include algebraic non- 
linear terms of power 2d + 1 in dimensions d |8|- Thus, 
an understanding of the joint effect of long-range interac- 
tions and of the agents' motion is desirable and relevant 
from ecological and statistical physics points of view. 

In this Rapid Communication, we aim to shed fur- 
ther light on the remarkable robustness of the LPPM 
scenario. To this end, we study an ecologically inspired 
stochastic lattice predator-prey model with next-nearest- 
neighbor interaction (NNN-LPPM), both in the presence 
and absence of a nearest- neighbor exchange process ('stir- 
ring'). We will demonstrate a subtle interplay between 
the correlations generated by the NNN interaction and 
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the stirring process. As a result, there is a regime where 
the NNN-LPPM phase diagram indeed follows the LPPM 
scenario outlined above, with a continuous predator ex- 
tinction transition in the DP universality class. On the 
other hand, we shall also see under which unexpected 
conditions a first-order phase transition can occur as a 
consequence of the competition between the short-range 
exchange and the NNN predator-prey interactions. 

To begin, we outline the main properties of the de- 
terministic LVM and then of the corresponding LPPM 
0,0,0- Consider two chemical species subject to the re- 
actions A — > (decay rate p > 0), B + — > B + B 
(branching rate a > 0), and A + B — > A + A (pre- 
dation rate A > 0). Neglecting any spatial variations 
and fluctuations of the concentrations a(x,t) and b(x,t) 
of 'predators' A and 'prey' B, one obtains the classi- 
cal LVM rate equations: a(t) = Xa(t)b(t) — pa(t) and 
b(t) = crb(t) — Xa(t)b(t). These deterministic equa- 
tions have as stationary states (a*,b*) — (0,0) (extinc- 
tion), (0, oo) (predators extinct, prey proliferation), and 
(a*,b*) — (a/X,p/X) (species coexistence). The un- 
stable fixed points (0, 0) and (0, oo) constitute absorb- 
ing states of the dynamics. The existence of a con- 
served first integral of the deterministic rate equations, 
K(t) = X[a(t) + b(t)] - crlna(t) - plnb(t) = constant, 
implies oscillatory kinetics around (a*,b*). 

Since this center singularity is unstable with respect 
to introducing model modifications ja, |jfl , the LVM rate 
equations are often rendered more 'realistic' by intro- 
ducing growth-limiting terms 0,0|- For the LVM, this 
amounts to replacing the rate equation for species B 
with b(t) = a b(t) [l - p- 1 b(t)] - X a(t) b(t) (p is the prey 
'carrying capacity'; growth-limiting terms for the preda- 
tors do not induce any qualitative changes.) The three 
fixed points are now shifted to (a*, b*) = (0,0) (extinc- 
tion), (0, p) (predators extinct, system saturated with 
prey), and (a*,b*) with a* = (1 — p/Xp) a/X, which is 
in the physical region (0 < a* < 1) if A > p/p, and 
b* = p/X. Linear stability analysis reveals (0,0) to be a 
saddle-point, whereas (0, p) is stable (node) if A < p/p 
(when a* < 0), and a saddle-point (stable in the b di- 
rection) otherwise. When A > p/p, the coexistence state 
(a*,b*) is stable; it is either a node or a focus, asso- 
ciated with spiral trajectories in the (a, b) phase plane 
0- Thus, at the rate equations level, A c = p/p is the 
critical predation rate. The global stability of (a*,b*) 
is established by the existence of a Lyapunov function 
C(a,b) = X{a*lna(t)+b*lnb{t)-a(t)-b(t)} Q. Many of 
these features re-emerge in stochastic LPPM with site re- 
striction. Indeed, Monte Carlo simulations fill f]~3| yield 
that as in mean-field theory the coexistence fixed point 
is either a node or a focus. In the latter case, amazingly 
rich spatio-temporal patterns of persistent predator-prey 
'pursuit and evasion' waves 0,0 emerge, inducing erratic 
correlated population density oscillations. In finite sys- 
tems, these quasi-periodic fluctuations appear on a global 
scale, but the amplitude of the density oscillations de- 
creases with system size [T5 | . A completely different pic- 



ture emerges when the active fixed point is a node, just 
above the predators' extinction threshold A c : Instead of 
the intricate front patterns, small predator 'clouds' effec- 
tively diffuse in a sea of prey ^3 . If the value of A is re- 
duced further (keeping the other rates fixed), at the criti- 
cal value A c the system reaches the absorbing state. This 
active-to-absorbing phase transition is found to be in the 
directed percolation (DP) universality class fl9l : this is 
also true for many LPPM variants [l3l . IT3. Il5| . These 
results can be understood from the master equation: For 
the above reactions one may derive an equivalent field 
theory action poj . which near A c can be mapped onto 
Reggeon field theory [la. known to describe the asymp- 
totic DP scaling laws [m IM l2l| . 

In most LPPM (see, e.g., Refs. [H El), the 'preda- 
tion' process subsumes nearest-neighbor interaction and 
the effects on both the prey and the predators in a sin- 
gle 'reaction'. More realistically, one should split this 
into two processes, and thereby introduce two indepen- 
dent time scales. This leads to the following stochas- 
tic reaction scheme that incorporates a three-site (NNN) 
process: (a) A predator reproduces in the vicinity of a 
prey (favorable environment) according to the triplet re- 
action A + (Z> + B^A + A + B [with rate S/z(z - 1); 
z = 2d is the coordination number of a d-dimensional 
hypercube] ; (b) a predator consumes a neighboring prey 
(rate r\/z), leaving an empty site, according to the bi- 
nary process A + B — ► + A; (c) we shall also allow 
for an efficient mixing process, through particle exchange 
with rate T>/z ('stirring') between two neighboring sites 
regardless of their content ^3- Besides these reactions, 
we still consider the processes B + — > B + B (rate 
a/z) and A — > (rate p). Assuming full site restriction, 
i.e. allowing at most one particle per site, the mean-field 
(MF) rate equations now read 

a(t) = S a(t) b(t) [1 - a{t) - b{t)\ - p a(t) (1) 
b(t) = a b(t) [1 - a(t) - b(t)] - f] a(t) b(t) . (2) 

These equations can be obtained from the master equa- 
tion of our NNN-LPPM upon factorizing the three-point 
correlators as products of the corresponding densities a, b. 
In contrast with the LVM, the nonlinear term in Eq. Q 
is cubic (NNN interaction); the site restriction appears 
through the growth-limiting factors 1 — a — b. Note that 
the mixing parameter T> does not enter the rate equations 
(but would appear in the equations for the three-point 
and higher correlation functions). Equations QJ, (0) ad- 
mit four fixed points, provided 6 > S c — Ap(a + f])/r). In 
addition to the previous absorbing states, (a*,b*) = (0, 0) 
and (0, 1), the two new nontrivial steady states (k = 1, 2) 
are given by 
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and b* k = |[1 + (-l) fe ^l - 5JS\. These active fixed 
points («! 21^12) correspond to two distinct predator- 
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FIG. 1: Flows in the phase plane from integrating Eqs. 
© for [i = a = 77 = 1, with 5 = 4 (left) and 8 = 9 (right). 
Left: (0, 1) is the only stable fixed point (node). Right: there 
is an additional stable (node) active fixed point (a*, &*) = 
(1/3,1/3); while (0,0) and (02,62) = (1/6,2/3) are unstable. 
The slopes of the separatrices at (02,62) are ~ 1.126 (dashed 
line) and » —1.568 -(see text). 
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FIG. 2: (Color online) Average stationary prey density b* vs 
S for 77 = a = 2/j = 2. Left: DP-like transitions on 256 2 , 50 3 
and 20 4 lattices when T> = 0. Right: Effect of the stirring 
on a 512 2 lattice. For T> = there is a continuous transition 
(center curve, black), while for T) = 10 (sufficient stirring) two 
stable branches emerge and there is a first-order transition; 
the top (red) branch corresponds to predator extinction, and 
the bottom (blue) one is associated with a coexistence phase. 



prey coexistence phases. From linear stability anal- 
ysis we infer that the absorbing state (0, 1) is al- 
ways a stable node: the associated Jacobian eigen- 
values read e+(0, 1) = — /u, [with eigenvector v + = 
({/1 — 0"}/{cr + i]}, 1)] and e_(0, 1) = —a [eigenvector 
d_ = (0,1)]. On the other hand, (0,0) is an unstable 
saddle-point, with eigenvalues e+(0, 0) = er [with unsta- 
ble eigendirection v + = (0,1)] and e_(0,l) = — // [sta- 
ble eigendirection V- = (1,0)]. Without loss of gen- 
erality, we just discuss the stability of the active fixed 
points |J3J when 77 = [i = a = 1. In this case, 8 C = 8 
and the eigenvalues of the Jacobian respectively read 

(k = 1,2): e±(a* k ,b* k ) = - 
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t v '22(-l) fc y/1 - 8/5 + 10 - 8/8. Thus, the active fixed 

point (a*,&!) is stable, $t(e±(al, &*)) < 0, while (og,^) 
is a saddle-point. More generally, for fixed /1, 17,77 there 
exists a value <5 S > <5 C such that (a*,b*) is a sta- 
ble node if S c < 8 < 8 S , and a stable focus, i.e., 
S(e±(a1,bl)) ^ 0, if 8 > 5 S . When rj = /j = a = 1, 
J, = (29 + ll\/7)/6 « 9.68388. Typical phase portraits 
as predicted by Eqs. Q, are illustrated in Fig. ^ 
The eigenvectors associated with (a%,b%) give the slope 

of the separatrices in its vicinity: 4(\f&— \f8 — 8)/ ^\/<5 — 

V^ - ^ ± y / 10<5 + 22^(5((5-8)-8) . It follows from this 

discussion that, at the mean- field level, the introduction 
of a triplet interaction changes the behavior of the system 
dramatically: For 6 > S c the system can reach the ab- 
sorbing state full of prey or alternatively a phase where 
the prey population, with stationary density b* < 1/2, 
coexists with the predators. Hence, the rate equations 
predict the possibility of a first- order phase transition. 

Motivated by these predictions, quite different from 
those of other LPPM, we have studied the properties of 
our NNN-LPPM through Monte Carlo simulations on pe- 
riodic hypercubic lattices. We have first considered the 
case of slow (T> « 0) and fast stirring and noticed the 
emergence of quite different behavior. In fact, for no (or 



slow) stirring, instead of a discontinuous phase transition, 
we have observed a continuous active-to-absorbing phase 
transition as for the LPPM in dimensions d — 2, 3 and 
even d = 4, see Fig. (left). (Of course, in dimensions 
d > 3 the model is biologically irrelevant: these cases 
have only been considered to assess the validity of the MF 
theory.) To ascertain the properties of the NNN-LPPM 
we have employed the dynamical Monte Carlo technique 
|l8j. Near the extinction threshold, one expects power 
law behavior for the survival probability P(t) ~ t~ s and 
the number of active sites N(t) ~ t 9 . By averaging over 
3 x 10 6 independent runs, performed on a 512 x 512 lat- 
tice, each with duration 10 5 Monte Carlo steps, for fixed 
rates 77 = a = 2/i = 2, and T> = we have estimated 
the critical point to be at S c rj 11.72 (the MF predic- 
tion is S c = 8), and measured 8' ~ 0.451 and 9 w 0.230, 
very close to the established two-dimensional DP expo- 
nents ■ A s illustrated in Fig. we have also deter- 
mined the order parameter critical exponent defined via 
a(t — > 00) — (8~8 c ) f3 as (3 w 0.584. We have checked that 
the exponent values are consistent with the DP univer- 
sality class for several choices of the rates 77, fi, a. Qual- 
itatively, the features of the NNN-LPPM remain similar 
in d = 3 and 4 (Fig.0 left): we observe continuous phase 
transitions (for different values of 8 C ) with (3 w 0.81 for 
d = 3 and f3 ~ 1.0 for d = 4 (upper critical dimension of 
DP) in agreement with DP values [l8j . 

In the absence of stirring, the phase diagram changes 
qualitatively when d > 5: even for T> = 0, one then ob- 
serves the first-order phase transition predicted by the 
MF approximation. The situation turns out to be com- 
pletely different when the stirring is sufficiently fast, as 
illustrated in Fig. (right): A first-order phase transi- 
tion occurs in low dimensions as well, and, depending on 
the initial condition with respect to the separatrices (see 
Fig. ^ right), the flows in the phase portrait end either 
at the absorbing fixed point (0,1), or reach a station- 
ary state where both predators and prey coexist (with 
b* < 1/2). This scenario, in the presence of sufficiently 
fast stirring, therefore recovers the mean-field behavior, 
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FIG. 3: (Color online) Average stationary density of preda- 
tors in the absence of stirring on a 512 x 512 lattice with 
r\ = a = 2/j, = 2 and T> = 0: existence of a DP-like phase 
transition at S c ~ 11.72 with exponent j3 ft* 0.584 (see inset). 

at least qualitatively. It is quite remarkable that the rate 
equations JTJ), 10 describe the NNN-LPPM already for 
mere nearest-neighbor (NN) exchanges at finite rates; one 
would rather expect the MF regime to emerge in the limit 
of infinitely fast exchange processes involving the swap of 
all particles (not restricted to NN partners) ^| • 

As illustrated in Fig.0] the intriguing properties of the 
NNN-LPPM with NN exchange process can be summa- 
rized as follow: (i) For vanishing mixing (T> small com- 
pared to the other rates), in dimensions 1 < d < 4 the 
system undergoes an active-to-absorbing state transition 
which belongs again to the DP universality class; only for 
d > 5, a first-order phase transition appears. Stochastic 
fluctuations clearly have a drastic effect here, invalidating 
the mean-field picture in dimensions d < 4. (ii) When one 
allows for random short-range particle mixing (T> > 0), 
the dynamics and the phase portrait flows change dra- 
matically [Fig. center], (iii) When the exchange pro- 
cesses become sufficiently fast (typically, when T> w S) 
a new fixed point associated with a coexistence phase 
is available (this holds even in d = 1), as demonstrated 
in Fig. 2] (right), and the system undergoes a first-order 
phase transition as predicted by the mean-field theory. 
As expected, when there is 'fast' stirring (T) much larger 



than the other rates) the MF predictions become very ac- 
curate. We have also checked that the NNN-LPPM sta- 
ble active fixed point is, in agreement with the MF anal- 
ysis and generic properties of the other LPPM 0, ^[ , 
either a node or a focus. When it is a focus, the coex- 
istence phase is again characterized by population oscil- 
lations originating in moving activity fronts but, as the 
system is more mixed, these 'rings' appear less prominent 
than in the LPPM with NN interactions ^(| . 

In this paper, we have first outlined the main prop- 
erties of the LPPM with nearest-neighbor interactions: 
namely, the existence of erratic oscillations and complex 
patterns deep in the coexistence phase and a directed 
percolation type phase transition. We have then further 
tested this scenario by considering a perhaps more re- 
alistic model variant with next-nearest-neighbor interac- 
tion. Upon in addition introducing a short-range stirring 
mechanism together with this longer-range interaction, 
an intriguing interplay emerges: When the NN exchange 
process is 'slow', the NNN reaction induces subtle corre- 
lations that completely invalidate the MF treatment and 
the system still undergoes a DP-type phase transition 
(for 1 < d < 4). In this regime, the generic LPPM sce- 
nario is thus fully confirmed. However, when the value 
of the mixing rate T> is raised, the simple NN exchange 
process 'washes out' the NNN correlations and the sys- 
tem reproduces the MF behavior, displaying a first-order 
phase transition. This is to be viewed in contrast with the 
standard LPPM, for which even fast diffusion of preda- 
tors and prey generally does not qualitatively affect its 
properties |14| . 
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